function [uIn] = hic(params)

  Nx = params.Nx;
  Ny = params.Ny;


  % ****************************
  % initial conditions
  % ****************************
  T  = params.Temp .* ones(Nx, Ny);
  ux = zeros(Nx, Ny) + 0.0001;
  uy = zeros(Nx, Ny) + 0.0001;
  P1 = zeros(Nx, Ny);
  P2 = zeros(Nx, Ny);

  % ********************************************************
  %     V O R T I C E S
  % ********************************************************

  % initial conditions 
  if params.hydroVortices
    vp = params.hydroVorticesParams;
    for i =1:size(vp,1)
      v = vp(i,:);
      vamp = v(4);
      vwidth = v(3);
      vx = v(1);
      vy = v(2);

      cx = vx.*params.Nx;
      cy = vy.*params.Ny;

      gauss_grid_width = vwidth .* params.Nx;
      sig2 = 2.*(gauss_grid_width).^2;

      for x = 1:params.Nx
        for y = 1:params.Nx
          for ddx = -1:1
            for ddy = -1:1
              %T(x, y) = T(x,y) + params.hydroPeakAmp .* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
              ux(x,y) = ux(x,y) -vamp.*(y+ddy*params.Ny-cy).*exp(-((x+ddx*params.Nx-cx)^2 + (y+ddy*params.Ny-cy)^2)/sig2);
              uy(x,y) = uy(x,y) +vamp.*(x+ddx*params.Nx-cx).*exp(-((x+ddx*params.Nx-cx)^2 + (y+ddy*params.Ny-cy)^2)/sig2);
            end
          end
        end
      end  % end of fors
    end
  end

  % ********************************************************
  %     E X P    P E A K
  % ********************************************************

  if params.hydroPeak
    cx = 0.5 * params.Nx + 0.5;
    cy = 0.5 * params.Ny + 0.5;
    cx2 = 0.5 * params.Nx + 0.5;
    cy2 = 0.6 * params.Ny + 0.5;
    gauss_grid_width = params.hydroPeakRelWidth .* params.Nx;

    sig2 = 2.*(gauss_grid_width).^2;
    for x = 1:params.Nx
      for y = 1:params.Nx
        T(x, y) = T(x,y) + params.hydroPeakAmp .* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
      end
    end
  end
  %ux = 0.99999999999999;

  [gx, gy] = ndgrid(params.Cx/params.Lx, params.Cy/params.Ly);
  gx = 2.*pi.*gx;
  gy = 2.*pi.*gy;


  % ********************************************************
  %     R A N D O M
  % ********************************************************

    rng(123);
  if params.hydroRandGradOn
    % smooth random gradients
    rng(123);
    N = 1000;
    %random_factor = 1/50; %correct for random interference
    for s=1:N
      phasex = rand(1)*2*pi;
      phasey = rand(1)*2*pi;
      ampx = rand(1);
      ampy = rand(1);
      freqx  = round(rand(1)*params.hydroRandGradFreq);
      freqy  = round(rand(1)*params.hydroRandGradFreq);
      ux = ux + 1/N.*ampx.*sin(freqx.*gx + phasex).*sin(freqy.*gy + phasey);

      phasex = rand(1)*2*pi;
      phasey = rand(1)*2*pi;
      ampx = rand(1)*params.hydroPureWaveAmp;
      freqx  = round(rand(1)*params.hydroRandGradFreq);
      freqy  = round(rand(1)*params.hydroRandGradFreq);
      uy = uy + 1/N.*ampx.*sin(freqx.*gx + phasex).*sin(freqy.*gy + phasey); %.*sin(freqy*gy + phasey);
    end
    maxx = max(max(abs(ux)));
    maxy = max(max(abs(ux)));
    disp(sprintf('max ux: %f  max uy: %f', maxx, maxy));

    if true
      normalization = 0.1;
      ux = ux / max([maxx maxy]) .* normalization;
      uy = uy / max([maxx maxy]) .* normalization;
    end
  end

    
  if true

    % turbulent flow
    ux = ux + 0.1 * sin(2*gy);%.*(1+0.2.*cos(2*gx));
    uy = uy + 0.005.*sin(2*gx+rand(1)*cos(3*gy));
    uy = uy + 0.003.*sin(1*gx+rand(1)*cos(2*gy));
    uy = uy + 0.0002.*sin(5*gx+rand(1)*cos(4*gx));
    amp = 0.1;
    P1  = P1 + amp*1*sin(4*gy+rand(1));
    P2  = P2 + amp*1*sin(3*gx+rand(1));
  end



  if params.hydroPureWave
    %ux = ux + params.hydroPureWaveAmp.*sin(params.hydroPureWaveFreq.*gy);
  end

  if false
    uxmax = max(max(abs(ux)));
    uymax = max(max(abs(uy)));
    umax = max([uxmax uymax]);
    ux = ux./umax * 0.5;
    uy = uy./umax * 0.5;
  end

  if false
    % copied from allan's file
    x = -pi + 2*pi*(0:Nx-1)'/Nx;
    y = -pi + 2*pi*(0:Ny-1)'/Ny;
    [X,Y]=ndgrid(x,y);

    % T = 0.*T + 1.0/50;
    ux = 0*X + 0.5 * sin(2*Y);
    uy = 0*X + 0.0001*sin(3*X+rand(1)*cos(7*Y));
    P1  = 0.000001*sin(3*X+rand(1));
    P2  = 0.000001*sin(2*Y+rand(1));



      cx = 0.5 .* params.Nx;
      cy = 0.3 .* params.Ny;
      sig2 = (0.05 .* params.Nx).^2 * 2;
      for x = 1:params.Nx
        for y = 1:params.Nx
          for ddx = -1:1
            for ddy = -1:1
              uy(x, y) = uy(x,y) + 0.05.* exp(-((x-cx)^2 + (y-cy)^2)/sig2);
              %ux(x,y) = ux(x,y) + rand*0.0001;;
              %uy(x,y) = uy(x,y) + rand*0.0001;;
            end
          end
        end
      end  % end of fors
  end



  % ****************************
  %  rescale temperature
  % ****************************

  % T = T .* (4 * pi / 3);



  % repackage initial conditions into a vector
  Nt = Nx.*Ny;
  uIn = zeros(Nt.*5, 1);
  uIn(0*Nt+1 : 1*Nt) = T(:);
  uIn(1*Nt+1 : 2*Nt) = ux(:);
  uIn(2*Nt+1 : 3*Nt) = uy(:);
  uIn(3*Nt+1 : 4*Nt) = P1(:);
  uIn(4*Nt+1 : 5*Nt) = P2(:);

end
